NUCLEAR SCIENCE AND TECHNIQUES 26, 050603 (2015) 


Monte Carlo fission matrix acceleration method with adaptive mesh* 


PAN Liu-Jun QAR), WANG Rui-Hong (£5 Z),' and JIANG Song ({L#) 
Institute of Applied Physics and Computational Mathematics of Beijing, 
No.2 Fenghao Donglu, Haidian District, Beijing 100094, China 
(Received April 29, 2014; accepted in revised form December 9, 2014; published online October 20, 2015) 


In large-scale, loosely coupled systems, Monte Carlo criticality calculation always suffers from slow fission 
source convergence resulting from the high dominance ratio. The fission matrix acceleration method has shown 
its potential to accelerate the convergence of the fission source in many numerical simulations. In practice, 
however, instability of this method may be caused by imbalanced precisions of elements of the fission matrix. 
Hence, an improved method, in which the space mesh used to compute the fission matrix is defined adaptively 
based on the fission bank in each cycle, is introduced. The proposed method ensures balanced precisions of 
elements of the fission matrix, so is more stable than the existing fission matrix method. 
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I. INTRODUCTION 


Monte Carlo (MC) method is widely used to compute the 
multiplication factor and fundamental mode eigenfunction 
of critical systems in reactor physics. The power iteration 
method is the most important technique for criticality calcu- 
lation in general Monte Carlo codes. However, this method 
always suffers from slow fission source convergence in large, 
loosely coupled systems. To achieve a sufficiently converged 
fission source, a large number of inactive cycles must be per- 
formed before active cycles start. This slow convergence 
may result in unacceptable computational cost, especially for 
whole core analysis of commercial reactors. To overcome 
this issue, various acceleration methods have been suggested 
in MC criticality calculations. The super-history method [1] 
and the Wielandt’s method [2] decrease the dominance ratio 
by modifying the power iteration scheme, and then reduce the 
necessary number of inactive cycles to reach the convergence 
of the fission source. However, they increase the time ex- 
pense of computation during each inactive cycle. The func- 
tional method [3] and the coarse-mesh finite difference ac- 
celerating Monte Carlo method [4] employ standard Monte 
Carlo techniques to estimate nonlinear functionals, which are 
used in low-order functional Monte Carlo (FMC) equations 
and coarse-mesh finite dfference (CMFD) equations to obtain 
the eigenvalues and discrete representations of the eigenfunc- 
tion. The modified power iteration method [5] can be used for 
convergence acceleration in addition to the capability of cal- 
culating the second eigenpair. There are still no examples or 
illustrations to demonstrate that the modified power iteration 
method will work for continuous energy problems. 

Compared with the methods mentioned above, the fission 
matrix acceleration (FM) method [6] is more intuitive. When 
the system is spatially discretized, the fission matrix can be 
tallied during the Monte Carlo simulation at each cycle. The 
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FM method adjusts the weights of fission neutrons according 
to the fundamental mode eigenvector of the fission matrix at 
each cycle. This method has shown its potential to acceler- 
ate the fission source in many numerical simulations. How- 
ever, some instability issues restrict the applications of this 
method. Actually, the FM method constructs a space mesh 
before the iteration process and uses this fixed mesh to tally 
the fission matrix in each cycle. The elements of the fission 
matrix may be computed with various precisions, as the el- 
ements correspond to zones with different source intensities 
are tallied with a different number of samples. The fission 
source intensities may become very small in some zones af- 
ter several cycles. This results in large statistical errors of the 
corresponding columns of the fission matrix. The imbalanced 
precisions of elements of the fission matrix may cause the FM 
method to be instable. This instability issue of the FM method 
was recently addressed in the semi-fixed-source FM method 
[7]. This method samples all elements of the fission ma- 
trix using an equal number of particles to improve efficiency. 
But the semi-fixed-source FM method increases the compu- 
tational cost, as it samples many more particles than the ex- 
isting method, especially when the real fission source of the 
system is severely imbalanced. In this paper, the space mesh 
used for computing the fission matrix is defined adaptively 
based on the fission bank in each cycle. This method ensures 
balanced precisions of elements of the fission matrix without 
increasing the number of sampled particles. 

It is worth noting that “accumulation” [8] can be adapted 
to address similar issues. This method improves the precision 
of the fission matrix, but it is not able to balance the precision 
of elements of the fission matrix. 


Il. FISSION MATRIX METHOD WITH ADAPTIVE MESH 
A. Fission matrix method 


The main purpose of the criticality calculation is to find the 
fundamental mode solution of the eigenvalue equation [7] 


Hs(r, E) = ks(r, E), (1) 
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where k is the eigenvalue, s(7, Æ) is the corresponding fission 
source distribution function, and 


Hs(F, B) = J dE! / di’ f(F, B! > F, E)s(7", E"), 
O V 
(2) 


in which f(r”, E” — r, E) is the fission kernel and represents 
the expected number of first generation fission neutrons pro- 
duced per unit volume at r per unit energy at E from a fission 
neutron born at 7” with energy E’. The fission neutrons are 
assumed to emit isotropically (this assumption is not neces- 
sary for Monte Carlo simulation, only for simplifying the for- 
mulation). 

As we know, the standard power method is used to solve 
problem (1) in most Monte Carlo codes. A guessed fission 
source distribution or a distribution obtained by other calcula- 
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The “cycle fission matrix” [9] can be constructed in each 
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As the fission source distribution evaluated from the fission 
matrix converges more rapidly than ordinary Monte Carlo 
simulations [7], the fission source convergence can be accel- 
erated by adjusting the weights of fission neutrons according 
to the fundamental mode eigenvector of H‘”) in each cycle. 
If the m-th fission neutron is produced in region 2, its weight, 


wi ) is corrected to [10] 


A~ Zi ym (6) 


where Ja is the fundamental eigenvector of H, and 


oF is the relative source intensity at zone i obtained from 
an ordinary MC simulation (normalized to the same number 


with So) In practice, the fission matrix acceleration is only 
used in inactive cycles. It is recommended that a few inac- 
tive cycles without the weight correction should be set be- 
fore moving into active cycles to eliminate the bias of region- 
averaged correction factors. 
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tions is used as an initial fission source distribution, and then 
the distribution is iteratively recalculated until it approaches 
to convergence. The cycles before obtaining a convergent 
fission source are called inactive cycles, while active cycles 
are the cycles after achieving a convergent fission source. In- 
active cycles are used to converge to a fission source and ac- 
tive cycles are used for tallying the final results. 

When the system is spatially discretized, the real fission 
source distribution satisfies the following equation: 


HS = kS, (3) 


where S indicates a fission source distribution vector, its el- 
ements represent the fraction of the fission source in each 
zone. k is the multiplication factor, and the fission matrix 
H is defined as follows: 


(4) 


m 


cycle as 


(5) 


n 


B. Adaptive mesh 


In Monte Carlo simulation, the number of histories sam- 
pled in each zone is decided by the source intensity. When a 
fixed mesh is used for computing the fission matrix during all 
cycles in the FM method, the fission source intensities may 
become very small in some zones after several cycles. There- 
fore, the number of histories sampled in those zones are less. 
Imbalanced precisions of elements of the fission matrix may 
occur and result in instability of the FM method. In this pa- 
per, to improve efficiency and overcome instability, adaptive 
mesh is proposed as follows: 


(1) A relatively fine mesh is defined, consisting of a relatively 
large number of zones. 


(2) At the beginning of each cycle, construct a coarse mesh 
by combining the zones of the fine mesh according to 
the source distribution, making the source intensities of 
different zones of the coarse mesh roughly equal. Each 
zone of coarse mesh may include several zones of the fine 
mesh. A different number of zones in the coarse meshes 
for different cycles is allowed. 


(3) During each cycle, the fission matrix is calculated based 
on the coarse mesh. Correspondingly, the order of the 
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fission matrixes may be different in different cycles. 


(4) During each cycle, count the number of the first gener- 
ation fission neutrons in each zone of the fine mesh and 
construct the fission bank for the next cycle. 


(5) After each cycle, adjust the weights of the neutrons in the 
fission bank according to the fundamental eigenvector of 
the fission matrix. 


Step (2) is the key step of this method. This step is easy to 
implement for a 1 dimensional model. For this purpose, the 
fine mesh contains N zones can be defined by 


a = £o < T1 <: < TN—1 < TN =D. (7) 


The source intensity (according to the fission source bank 
given by the previous cycle) of zone J; = [x;_1, xi] is As;. If 
we want to construct a coarse mesh containing about K zones 
(K should be much smaller than N ), the source intensity of 
each zone of the coarse mesh should be about 


N 
> As; 
i=1 
t= =... 8 
4 (8) 
We first select kı satisfying the following equation: 
kı k 
> ASit = min > ASit , (9) 
i=1 t=1 
kı 
let |]J J; be the first zone of the coarse mesh, then find kə 
i=1 
satisfying 
k2 k 
N ASi -t = min N ASit è, aO) 
i=kı i=kı 
k2 


let U JI, bethe second zone of the coarse mesh, and so on. 
i=kı+1 

After this, 1f the intensity of the last zone of the coarse mesh 
is less than a settled value (for example, t/2 ), the last two 
zones of the coarse mesh are combined. This method is not 
able to make the intensity of each zone of coarse mesh equal 
to t exactly, so the number of zones in the coarse meshes may 
be a little different than K . 


HI. NUMERICAL RESULTS 


We consider a one dimensional problem of large hetero- 
geneous fissile regions surrounded by a thin reflector. The 
neutrons are assumed to be scattered isotropically, and the 
physical data is given in Table 1. 

Monte Carlo simulations were performed with 500 000 his- 
tories at each cycle. The initial fission source distribution is 
spatially flat. The FM method is performed with the fissile 
regions discretized to 10 zones. The 30cm fissile region 
is equally discretized to 300 zones for the fine mesh. The 
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TABLE 1. Data for test problem 


Region Location (cm) SJ, (cm1) a (mT) vY. (cm?) 
1 O<x4<5 1.0 0.86 0.00 

2 5<x<10 1.0 0.86 0.15 

3 10 < x < 30 1.0 0.86 0.10 

4 30 < x < 35 1.0 0.86 0.149 

5 35 <x < 40 1.0 0.86 0.00 


amounts of zones of coarse mesh constructed for the fission 
matrix acceleration method with the adaptive mesh (ADFM) 
method are also about 10, but may be a little different at some 
cycles. 
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Fig. 1. The numbers of fission neutrons of each zone corresponding 
to different cycles when fixed mesh is used. 


Figure | shows the number of fission neutrons in each zone 
corresponding to the initial source, cycle 5 and cycle 20, re- 
spectively, when the fixed mesh is used. From this figure, it 
can be seen that the number of fission neutrons becomes very 
imbalanced after several cycles. This would result in very im- 
balanced precisions of elements of the fission matrix, because 
only the histories started from zone 7 contribute to the tallies 
of the elements H;; of the fission matrix. 
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Fig. 2. The coarse meshes for ADFM method corresponding to 
different cycles. 


Figure 2 shows the coarse meshes of the ADFM method 
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corresponding to different cycles. The mesh of the FM 
method is the same as the mesh of the ADFM method at the 
first cycle. Notice that the meshes are obviously different at 
different cycles. 

100 cycles are done with the FM method and the ADFM 
separately, the Shannon entropy [11] evolution is shown in 
Fig. 3 (the fine mesh is used for calculating the Shannon en- 
tropy). The average of the Shannon entropy over cycles 3001 
to 5000 of the standard Monte Carlo power iteration (PI) is 
taken as the reference. 
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Fig. 3. (Color online) Shannon entropies of first 100 cycles of 


different methods for one dimensional slab problem. 


In Fig. 3, it can be seen that both the Shannon entropies 
of the FM method and the ADFM method cross the refer- 
ence more rapidly than the standard PI and the Shannon en- 
tropies of the FM oscillate much more wildly than the ADFM 
method. 

In some situations, it may need many cycles to eliminate 
this oscillation after terminating the weight correction when 
the FM method is used. For example, the Shannon evolution 
when the weight correction is terminated at cycle 48 is shown 
in Fig. 4. 


standard PI 
/ 
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Shannon entropy 
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Fig. 4. (Color online) Shannon entropy evolution when the weight 
correction is terminated at cycle 48. 
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Table 2 shows the cycle numbers when the Shannon en- 
tropy crosses the reference for the first time after terminating 
the weight correction by FM and ADFM at different cycles 
(by standard PI, the Shannon entropy crosses the reference 
for the first time at cycle 851). It can be seen from this table 
that the ADFM method is more efficient and more stable than 
the existing FM method. 

The eigenfunction estimated during cycles 50 to 100 by the 
FM method, the ADFM method, and standard power iteration 
are shown in Fig. 5, in which the eigenfunction estimated dur- 
ing cycles 3001 to 3050 by the standard PI is taken as a refer- 
ence. 


FM at first 47 cycles 


© 
= standard PI 


ADFM at first 47 cycles 


Fig. 5. (Color online) Eigenfunction estimated during cycles 50 to 
100 of one dimensional slab problem. 


IV. DISCUSSION 


The convergence of source distribution can be diagnosed 
by checking whether the Shannon entropy crosses the refer- 
ence [11]. It can be seen that both the Shannon entropies 
of the FM method and the ADFM method cross the reference 
rapidly from the numerical results (Fig. 3 and Fig. 4). It seems 
that the convergence rate of the ADFM method is the same as 
the FM method. But it can also be seen that the FM method 
can be too oscillatory to use in this model. To eliminate this 
oscillation, extra cycles without weight correction should be 
set before moving into active cycles. Actually, when the FM 
method is used, terminating the weight correction at different 
cycles may lead to different effects. In some cases, it requires 
many cycles to make the Shannon entropy cross the reference 
again after terminating the weight correction. The amount of 
extra cycles may even be comparable to the amount of cycles 
of standard PI to make the Shannon entropy cross the refer- 
ence. But when ADFM is used, it always needs much fewer 
cycles to make the Shannon entropy cross the reference again 
after terminating the weight correction than standard PI (Ta- 
ble 2). Therefore, the ADFM method is more efficient and 
stable than the FM method. 

Compared to the FM method, the extra time required for 
the ADFM method can be neglected in the tested 1D model. 
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TABLE 2. Cycle numbers when the Shannon entropy crosses the reference at first time after terminating the weight correction 


The cycle terminating the weight correction 
The first cycle when the Shannon entropy crosses the FM 
reference after terminating the weight correction ADFM 


The same situation can be expected for 2D or 3D problems. 
Actually, the main extra cost of the ADFM method is to cal- 
culate the intensity of the fission source in each zone of the 
fine mesh. When ADFM is used, regular fine mesh can be 
defined (the union of the zones of fine mesh must contain the 
fissile regions but not need to be the same as the fissile re- 
gions). Then it is easy to find where a fission neutron in the 
fission source bank is located, so it is easy to count the num- 
ber of fission neutrons (which is the intensity of the fission 
source) in each zone of the fine mesh. 

The FM method has been implemented into the MCNP 
code to simulate 2D and 3D continuous-energy prob- 
lems [12]. The only difference between our ADFM method 
and the FM method lies in how to construct the space mesh. 
Although the ADFM method is only tested on a multi-group 
problem for simplicity, there is no difficulty in applying the 
ADFM method to continuous-energy problems. 

There is also no difficulty in applying the ADFM method 
to 2D and 3D problems. But for complex 2D or 3D systems, 
the coarse meshes must be constructed elegantly to make the 
ADFM method efficient. The reason is as follows. The 
weight corrections of the FM method and the ADFM method 
can only accelerate the balance of the fission source between 
different zones; they are not able to accelerate the balance of 
the local fission source in each zone. In the ADFM method, 
as the source intensities in each zone of the coarse mesh are 
roughly equal, local balance of the fission source in large 
zones are not important for their low source densities. For 1D 
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problems, the balances of local fission source in small zones 
are easy to reach. But for 2D or 3D problems, the small area 
or volume of a zone does not certainly lead to a quick bal- 
ance of the local fission source, because even 1f the area or 
volume of a zone is small, it is possible that the scale of one 
direction of that zone is large (for example, a long and narrow 
rectangle). To make the ADFM method efficient, the zones 
of coarse mesh must meet one extra requirement: the scale of 
one direction can not be much larger than others. 


V. CONCLUSION 


The fission matrix acceleration method has shown its po- 
tential to accelerate the convergence of the fission source in 
many numerical simulations. In practice, however, instability 
of this method may be caused by the imbalanced precision of 
elements of the fission matrix. In this paper, the space mesh 
used to compute the fission matrix is defined adaptively based 
on the fission bank in each cycle. This method ensures a bal- 
anced precision of elements of the fission matrix, so is more 
stable. Although the ADFM method is only tested on one 
dimensional multi-group problem for simplicity, there is no 
difficulty in using the ADFM method in two dimensional or 
three dimensional continuous-energy problems, on the condi- 
tion of properly constructed adaptive mesh. 
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